Settlement Persistence Project, SBOM Script #5:
Topography, Environment, Agriculture and Subsistence
Rudolf Cesaretti
Last run on 2022-07-09
- 1 Setup
- 2 Catchment Area
Topographic/Environmental Metrics
- 2.0.1 Potential Productivity
- 2.0.2 Subsistence/Agricultural Zones
- 2.0.3 Terrain Ruggedness
- 2.0.4 Hydraulic Agriculture Potential
- 2.0.5 Intensification Costs
- 2.0.6 Agricultural Potential
- 2.0.7 Population Pressure
- 2.0.8 Erosion Potential
- 2.1 Calculate Metrics for Catchments/Periods
- 2.2 Reorganize and Export Raster Data
- 3 Recombining and Reorganizing the Data
- 4 Export Data for Script #6
- References
I do four things in this R markdown document: Topographic/environmental metrics
- Calculate catchment area variables, including
- Topographic/environmental metrics
- Reorganize the data and export for Script #6
1 Setup
All of the data and scripts are downloadable from the new ASU SettlementPersist2022 github repository, which can be downloaded locally as a .zip folder or cloned to your own account.
Either way, once you have done so, you will need to modify the working directory (setwd(“C:/…)”) path and “dir” variables in the code chunk below to match the repository location on your computer.
wd <- list()
# SET YOUR LOCAL DIRECTORY LOCATION HERE:
wd$dir <- "C:/Users/rcesaret/Dropbox (ASU)/ASUSettlementPersist2022/"
# wd$dir <- 'C:/Users/TJ McMote/Dropbox (ASU)/ASUSettlementPersist2022'
wd$analysis <- paste0(wd$dir, "analysis/")
wd$data_r <- paste0(wd$dir, "data-raw/")
wd$data_p <- paste0(wd$dir, "data-processed/")
wd$data_f <- paste0(wd$dir, "data-final-outputs/")
wd$figs <- paste0(wd$dir, "figures/")
wd$funcs <- paste0(wd$dir, "functions/")1.1 Load R Packages and Custom Functions
# Package names
packages <- c("rgdal", "rgeos", "sp", "sf", "GISTools", "raster", "Matrix",
"gdistance", "lwgeom", "tidyverse", "tidyr", "stars", "dismo", "purrr",
"spatialEco", "whitebox", "classInt", "ggnewscale") #, 'data.table', 'zoo', 'era', 'JOPS', 'mgcv','igraph', 'ggrepel','ggridges', 'movecost', 'datplot', 'scales',
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# load packages
invisible(lapply(packages, library, character.only = TRUE))
rm(packages, installed_packages)
# Read in custom R functions located in the wd$funcs directory folder
FUNCS <- list("splitByAttributes.R", "RescaleSpatRast.R")
invisible(lapply(FUNCS, function(x) source(paste0(wd$funcs, x))))
rm(FUNCS)1.2 Import Data and Reorganize Data
Data we are importing:
- the AggSite polygon and catchment polygon data from script 4 (SPDFs)
- A simple polygon calculated in QGIS that specifies a hard outer border for the catchment areas of sites (constructed for sensitivity to survey borders and sites not included in the SBOM sample)
- Cost-distance matrices from script #3
- A raster hillshade basemap for the SBOM which includes the lakes
# Agg Site and catchment polygon data
All_AggPoly <- readOGR(paste0(wd$data_p, "SBOM_AggSitePoly4.gpkg"))## OGR data source with driver: GPKG
## Source: "C:\Users\rcesaret\Dropbox (ASU)\ASUSettlementPersist2022\data-processed\SBOM_AggSitePoly4.gpkg", layer: "SBOM_AggSitePoly4"
## with 1417 features
## It has 260 fieldsAll_CatchPoly <- readOGR(paste0(wd$data_p, "SBOM_CatchPoly4.gpkg"))## OGR data source with driver: GPKG
## Source: "C:\Users\rcesaret\Dropbox (ASU)\ASUSettlementPersist2022\data-processed\SBOM_CatchPoly4.gpkg", layer: "SBOM_CatchPoly4"
## with 1417 features
## It has 260 fields# reorder spatial points dataframe by period
All_AggPoly <- All_AggPoly[order(All_AggPoly$PeriodNum), ]
All_CatchPoly <- All_CatchPoly[order(All_CatchPoly$PeriodNum), ]
# Split polygons by Phase, saved as list of SPDF
Poly_List <- splitByAttributes(spdata = All_AggPoly, attr = "Period", suffix = "_SitePoly")
Catch_List <- splitByAttributes(spdata = All_CatchPoly, attr = "Period", suffix = "_CatchPoly")
# Catchment boundary limit polygon
CatchLims <- readOGR(paste0(wd$data_r, "CatchLims.gpkg"))## OGR data source with driver: GPKG
## Source: "C:\Users\rcesaret\Dropbox (ASU)\ASUSettlementPersist2022\data-raw\CatchLims.gpkg", layer: "CatchLims"
## with 1 features
## It has 89 fieldsCatchLims.sf = st_as_sf(CatchLims) #for ggplot
## Hillshade Basemap Raster with lake
HillshadeLake <- raster(paste0(wd$data_r, "HillshadeLake.tif"))
HillshadeLake <- rast(HillshadeLake, crs = 26914)
Hillshade.s <- st_as_stars(HillshadeLake) #for ggplot basemap
# Cost-distance matrices
ord <- c("EF", "EF_MF", "MF", "MF_LF", "LF", "LF_TF", "TF", "TF_CL", "CL", "CL_ET",
"ET", "ET_LTAzI", "LTAzI", "LTAzI_EA", "EA", "EA_LA", "LA")
temp = list.files(path = paste0(wd$data_p, "CDMatrices/"), full.names = TRUE)
nam.distmat = list.files(path = paste0(wd$data_p, "CDMatrices/"), full.names = F)
nam.distmat = gsub("_cdmat.csv", "", nam.distmat)
CD.mats = lapply(temp, read.csv, header = TRUE, row.names = 1)
CD.mats = lapply(CD.mats, as.matrix)
names(CD.mats) <- nam.distmat
CD.mats <- CD.mats[ord]
names(CD.mats) <- paste0(ord, "_cdmat")2 Catchment Area Topographic/Environmental Metrics
The environmental characteristics of site catchment areas may be a causally important factor in long-term settlement persistence. Two major environmental processes that might contribute to settlement persistence include subsistence/agricultural productivity and landscape degradation. We have therefore calculated a number of targeted environmental proxy metrics for these theoretical dimensions, outlined below. All of these metrics use only online environmental data and superficial contextual info about the study region
# create temp directory for rasters we will delete later we will stack and
# save these together at the end of this chunk
dir.create(paste0(wd$data_p, "temp_rasts"))
# sf_use_s2(FALSE)#Spherical geometry (s2) switched off
## Import 90 meter DEM
DEM <- raster(paste0(wd$data_r, "SBOM_DEM_90m.tif"))
DEM <- rast(DEM)
DEM <- project(DEM, "epsg:26914")
## Import SBOM Lake rasters calculated in GRASS GIS
LakesBinaryInv <- rast(paste0(wd$data_r, "LakesBinaryInv.tif"))
LakesBinaryInv <- project(LakesBinaryInv, DEM)
LakesBinaryInvNA <- LakesBinaryInv
LakesBinaryInvNA[LakesBinaryInvNA == 0] <- NA
LakesBinaryNA <- LakesBinaryInv + 1
LakesBinaryNA[LakesBinaryNA == 2] <- NA
LakesBinaryNA.s <- st_as_stars(LakesBinaryNA)
lake.s = st_as_sf(LakesBinaryNA.s, merge = TRUE)
SaltLakeBinaryInv <- rast(paste0(wd$data_r, "SaltLakeBinaryInv.tif"))
SaltLakeBinaryInv <- project(SaltLakeBinaryInv, DEM)
SaltLakeBinaryInvNA <- SaltLakeBinaryInv
SaltLakeBinaryInvNA[SaltLakeBinaryInvNA == 0] <- NA
FreshLakeDepth <- rast(paste0(wd$data_r, "FreshLakeDepth.tif"))
FreshLakeDepth <- project(FreshLakeDepth, DEM)2.0.1 Potential Productivity
As the basis and constraint of all ecosystem services, Net Primary Productivity (NPP) here serves as a basal proxy for potential subsistence/agricultural productivity without any landesque capital intensification. NPP is the amount of plant biomass per unit area minus respiratory costs (Gross Primary Productivity = GPP = total biomass per unit area).
Here, NPP is calculated for the SBOM using Nick Gauthier’s script “Climatic Potential for Agricultural Productivity” (npp_plots.Rmd) located in Matt Peeples’ Settlement Persistence GitHub repository via proxy calculation from precipitation and temperature. These raster inputs for Nick’s script (cropped to the SBOM study area) are included in the Script #3 directory folder.
TEXT FROM NICK GAUTHIER: The Miami model of Helmut Lieth characterizes climatic potential net primary production (NPP) as a function of two limiting factors: mean annual temperature and total annual precipitation. It uses the following empirical functions to relate each variable to potential NPP, the lowest of which (i.e. the limiting factor) is taken as the final NPP estimate:
\[NPP_{precip} = 3000(1 - e^{-0.000664
P})\] \[NPP_{temp} = \frac{3000}{1 +
e^{1.315-0.119T}}\] We’ll apply this model to 30-year
temperature and precipitation averages from the CHELSA dataset, a set of
statistically downscaled climate data that corrects for topoclimatic
influences like orographic precipitation and temperature inversions.In
order to run this code, you will need to download the large CHELSA
dataset files “CHELSA_bio10_01.tif” and “CHELSA_bio10_12.tif” from the
CHELSA server:
The resulting map gives us climatic potential NPP in grams of dry matter per square meter per year, with a theoretical maximum value of 3,000 \(g/m^2\). This can be easily converted to expected crop yields using empirical coefficients for different crop types.
The metrics calculated here are:
- NPP.tot = the sum of NPP in a catchment area. This is extensive metric is intended to capture the total scale/magnitude of resources available in a catchment area, similar to (a proxy for) carrying capacity without any landesque intensification.
- NPP.avg = Total NPP divided by the catchment area in hectares. This is intended to be an intensive proxy metric of a catchment’s subsistence productivity without any landesque intensification.
######### NPP #########
## Calculate net primary productivity (NPP)
temperature <- read_stars(paste0(wd$data_r,'CHELSA_bio10_01_SBOM.tif'))
precipitation <- read_stars(paste0(wd$data_r,'CHELSA_bio10_12_SBOM.tif'))
npp2 <- c(temperature, precipitation)
names(npp2) <- c("temperature", "precipitation")
npp <- npp2 %>% mutate(temperature = temperature / 10, # temperature is in degrees C * 10
### the miami model
npp_prec = 3000 * (1 - exp(-0.000664 * precipitation)),
npp_temp = 3000 / (1 + exp(1.315 - 0.119 * temperature)),
npp = pmin(npp_prec, npp_temp)) %>% select(npp)
rm(temperature,precipitation,npp_prec,npp_temp)
lake.s <- st_transform(lake.s, crs = 26914)
CatchLims.sf <- st_transform(CatchLims.sf, crs = 26914)
NPP <- st_transform(npp, crs = 26914)
# convert from stars to raster
NPP <- as(npp, "Raster")
NPP <- rast(NPP)
NPP <- project(NPP, "epsg:26914")
#Downscale to DEM resolution using DEM
NPP <- terra::resample(NPP, DEM, method="bilinear")
stz = st_as_stars(NPP)
plt = ggplot() +geom_stars(data = stz)+
scale_fill_viridis_c(na.value = NA, name = expression(g/m^2)) +
geom_sf(data = lake.s, fill="turquoise1", color="turquoise1", alpha = 0.2) +
geom_sf(data = CatchLims.sf, color="white", size=1, alpha = 0) +
ggtitle('Potential Unintensified Subsistence Productivity', 'Proxied by Maximum Potential Net Primary Production (NPP)') +
theme_bw() +
coord_sf(datum = sf::st_crs(26914))+
labs(x = 'Easting', y = 'Northing')
ggsave("NPPRast.png", plot = plt, device = "png", path = wd$figs, scale = 1, width = 6, height = 5, units = "in", dpi = 1000)
#NPP <- NPP * SaltLakeBinaryInvNA
NPP.rs <- RescaleSpatRast(NPP) # rescale 0-1
#NPP.rs <- NPP.rs+0.5# rescale 0.5-1.5
plt
2.0.2 Subsistence/Agricultural Zones
Sanders, Parsons & Santley (1979) recognize five major environmental zones in the Southern Basin of Mexico that strongly constrain agricultural practices in the region. Due to the impact of elevation on rainfall, frosts, topographic suitability, vegetation and soils, these environmental zones can be approximated from elevation alone:
- Lakebed Inundation Zone (<2240 masl) = 0
- Alluvial Plain (2240-2275 masl) = 1
- Lower Piedmont (2275-2350 masl) = 2
- Middle Piedmont (2350-2500 masl) = 3
- Upper Piedmont (2500-2700 masl) = 4
- Sierra (2700-6000 masl) = 5
Here, we calculate these zones from a 90m DEM.
The metrics calculated here are: * EZ.avg = Mean of ordinal numeric codes (see above) * EZ.sd = StDev of ordinal numeric codes (see above)
######### EnvZones #########
## Calculate BOM Environmental Zones from elevation (DEM)
m <- c(0, 2240, 0, # Inundation Zone / Lakebed Alluvium
2240, 2275, 1, # Lakeshore Alluvial Plain
2275, 2350, 2, # Lower Piedmont
2350, 2500, 3, # Middle Piedmont
2500, 2700, 4, # Upper Piedmont
2700, 6000, 5) # Sierra
EZm <- matrix(m, ncol=3, byrow=TRUE) #make into a matrix
EZ <- terra::classify(DEM, EZm, include.lowest=TRUE) # reclassify groups of values
stz = st_as_stars(EZ)
plt=ggplot() +geom_stars(data = stz)+
scale_fill_gradientn(colors=terrain.colors(6), values=c(0,0.2,0.4,0.6,0.8,1.2), na.value = NA) +
geom_sf(data = lake.s, fill="turquoise1", color="turquoise1", alpha = 0.4) +
geom_sf(data = CatchLims.sf, color="black", size=1, alpha = 0) +
ggtitle('Elevation-Conditioned Ecological Zones in the SBOM', 'Determinants of Subsistence/Agricultural Practices') +
coord_sf(datum = sf::st_crs(26914))+
theme_bw() +
labs(x = 'Easting', y = 'Northing')
ggsave("EnvZonesRast.png", plot = plt, device = "png", path = wd$figs, scale = 1, width = 6, height = 5, units = "in", dpi = 1000)
plt
2.0.3 Terrain Ruggedness
A major constraint of agriculture in any region is the suitability of local topography to arable cultivation. More steeply-sloping topography requires heavy labor inputs in the form of terracing to create arable surfaces and constrain erosion/runoff. By contrast, flat alluvial areas generally lack these de facto labor constraints, and are often better watered and blanketed in fertile and easily tilled sediment.
To capture this dimension, we use rescaled Terrain Ruggedness Index (TRI), which is a local (moving window), multi-directional, dimensionless index of local slope. More specifically, TRI is the mean of the absolute differences in elevation among a DEM raster cell and its 8 surrounding cells. Here, TRI is calculated using the “terrain” function in the “terra” R package, rescaled to the range [0,1]. Here, we calculate TRI in R using the “terra” package (Hijmans, 2022).
The metrics calculated here are:
- TRI.tot = Total (catchment sum) TRI = extensive metric intended to capture the total scale/magnitude of land that may need to be terraced
- TRI.avg = Mean TRI = intensive metric intended to capture the mean catchment area need for terracing
######### TRI #########
TRI <- terra::terrain(DEM, v = "TRI") # calc Terrain Ruggedness Index
# TRI <- TRI * SaltLakeBinaryInvNA
TRI.01 <- RescaleSpatRast(TRI) # rescale 0-1
stz = st_as_stars(TRI)
plt = ggplot() + geom_stars(data = stz) + scale_fill_viridis_c(option = "inferno",
na.value = NA, name = "TRI") + geom_sf(data = lake.s, fill = "turquoise1",
color = "turquoise1", alpha = 0.3) + geom_sf(data = CatchLims.sf, color = "white",
size = 1, alpha = 0) + ggtitle("Topographic Need for Terracing in the SBOM",
"Proxied by the Terrain Ruggedness Index (TRI)") + coord_sf(datum = sf::st_crs(26914)) +
theme_bw() + labs(x = "Easting", y = "Northing")
ggsave("TRIRast.png", plot = plt, device = "png", path = wd$figs, scale = 1,
width = 6, height = 5, units = "in", dpi = 1000)
plt
2.0.4 Hydraulic Agriculture Potential
The ability to increase agricultural output via irrigation and wetland agriculture landesque capital intensification (i.e. hydraulic engineering) is a critical determinant of socioeconomic performance across agrarian economies.
The suitability of an area to irrigation is controlled by the magnitude of locally available fluvial discharge, which is strongly conditioned by the topography of stream networks and drainage basins. All else being equal, the magnitude of fluvial discharge available to cultivators from the local drainage network is proportional to upstream drainage area (Anderson & Anderson, 2010). e proxy for this topographic potential drainage density is the Topographic Wetness Index, TWI. TWI measures the potential hydrological response of basins based on slope and the size of upslope contributing areas. This indicates which parts of the landscape are more likely to have overland flow, streams and saturation (Beven & Kirkby, 1979). Here, we calculate TWI in R using the “whitebox” package (Lindsay, 2016).
The suitability of an area to wetland agriculture is contexually dependent by region. In the SBOM, this was the area of the large, shallow freshwater lakes Chalco and Xochimilco, on which chinampas were built. The approx area of these lakes was modelled in GIS using published data and exported as a binary raster [0,1].
- IrrigPot.tot = Total Irrigation Potential = sum TWI
- IrrigPot.avg = Mean Irrigation Potential = sum TWI divided by catchment area
- WetAgPot.tot = Total Wetland Agriculture Potential = sum freshwater lake raster in catchment
- WetAgPot.avg = Mean Wetland Agriculture = sum freshwater lake divided by catchment area
######### TWI #########
wbt_breach_depressions(dem = paste0(wd$data_r, "SBOM_DEM_90m.tif"), output = paste0(wd$data_p,
"temp_rasts/DEM_breach.tif"))
DEM_breach <- rast(paste0(wd$data_p, "temp_rasts/DEM_breach.tif"))
DEM_breach <- project(DEM_breach, DEM)
writeRaster(DEM_breach, paste0(wd$data_p, "temp_rasts/DEM_breach.tif"), overwrite = TRUE)
# DEM2 <- as(DEM, 'Raster')
slope <- terra::terrain(DEM, v = "slope", unit = "degrees") #Slope
slope <- project(slope, DEM)
writeRaster(slope, paste0(wd$data_p, "temp_rasts/Slope.tif"), overwrite = TRUE)
wbt_d8_flow_accumulation(paste0(wd$data_p, "temp_rasts/DEM_breach.tif"), output = paste0(wd$data_p,
"temp_rasts/Accum.tif"), out_type = "specific contributing area")
Accum <- rast(paste0(wd$data_p, "temp_rasts/Accum.tif"))
Accum <- project(Accum, DEM)
Accum[is.nan(Accum)] <- NA
writeRaster(Accum, paste0(wd$data_p, "temp_rasts/Accum.tif"), overwrite = TRUE)
wbt_wetness_index(sca = paste0(wd$data_p, "temp_rasts/Accum.tif"), slope = paste0(wd$data_p,
"temp_rasts/Slope.tif"), output = paste0(wd$data_p, "temp_rasts/TWI.tif"),
verbose_mode = FALSE)
TWI <- rast(paste0(wd$data_p, "temp_rasts/TWI.tif"))
TWI <- project(TWI, DEM)
writeRaster(TWI, paste0(wd$data_p, "temp_rasts/TWI.tif"), overwrite = TRUE)
# TWI <- TWI * SaltLakeBinaryInvNA
TWI.01 <- RescaleSpatRast(TWI) # rescale 0-1
TWI.12 <- TWI.01 + 1
######### WetAG #########
WetAgPot <- rast(paste0(wd$data_r, "FreshLakeBinary.tif")) #binary 0,1 raster of lake
WetAgPot <- project(WetAgPot, DEM)
WetAgPot2 <- WetAgPot * 2
WetAgInv <- rast(paste0(wd$data_r, "FreshLakeBinaryInv.tif")) #binary wetAg area=0, other is 1
WetAgInv <- project(WetAgInv, DEM)
stz = st_as_stars(TWI)
plt = ggplot() + geom_stars(data = stz) + scale_fill_viridis_c(option = "mako",
na.value = NA, name = "TWI") + geom_sf(data = lake.s, fill = "turquoise1",
color = "turquoise1", alpha = 0.1, size = 0) + geom_sf(data = CatchLims.sf,
color = "grey", size = 0.5, alpha = 0) + ggtitle("Irrigation Agriculture Potential in the SBOM",
"Proxied by the Topographic Wetness Index (TWI)") + coord_sf(datum = sf::st_crs(26914)) +
theme_bw() + labs(x = "Easting", y = "Northing")
ggsave("TWIRast.png", plot = plt, device = "png", path = wd$figs, scale = 1,
width = 6, height = 5, units = "in", dpi = 1000)
plt
2.0.5 Intensification Costs
The costs of agricultural intensification are another possible factor in settlement persistence. The cost of terracing should be proportional to TRI, as noted above.
The cost of irrigation should be proportional erosive potential of fluvial discharge. Here, the crucial topographic variables are channel slope, \(S\), and upstream drainage area, \(A\) (assumed to be proportional to discharge). Assuming that tectonic uplift is negligible, the rate of erosion of any point in a drainage network is thus given by the stream power equation \[ \frac{\partial Z}{\partial t}=k_{sp}A_{sp}^{m}S_{sp}^{n} \] where \(Z\) is the depth of erosion, \(t\) is time, \(k\) is the scalar coefficient of erodibility, and \(m\) and \(n\) are positive dimensionless parameters that govern the relative importance of discharge to local slope for bed incision via shear stress (Anderson & Anderson, 2010; Whipple, 2004; Whipple & Tucker, 1999). This equation is normally used to model long-term channel incision and knickpoint migration in detachment-limited rivers, but we can also it’s general principles to measure topographic controls on the erosive potential of runoff-driven flood discharge in the seasonal gullys, torrents, and arroyos of the SBOM. The Stream Power Index (SPI) provides a proxy for this topographic potential for stream erosion (Moore et al., 1993), which we calculate in R using the “whitebox” package (Lindsay, 2016).
Finally, the cost of wetland agriculture should be proportional to freshwater lake depth. As such, we can proxy it by taking the estimated lake depth (elevation of lake area minus approx lake level of 2240 masl).
The relative_costs of these strategies are assumed to be 1:1:2 Terrace:Irrig:Wetland. As such, the rasters of these proxies will be rescaled to [0,1] for Terracing and Irrigation, and [1,2] for wetland.
The metrics calculated here are:
- IntnsCost.tot = Total catchment intensification costs = sum(TRI[0,1]) + sum(SPI[0,1]) + sum(LakeDepth[1,2])
- IntnsCost.avg = Average catchment intensification costs = sum(TRI[0,1]) + sum(SPI[0,1]) + sum(LakeDepth[1,2])/Area
- IntnsCostNW.tot = Total catchment intensification costs excluding wetland agriculture = sum(TRI[0,1]) + sum(SPI[0,1])
- IntnsCostNW.avg = Average catchment intensification costs = sum(TRI[0,1]) + sum(SPI[0,1])/Area
######### SPI #########
wbt_stream_power_index(sca = paste0(wd$data_p, "temp_rasts/Accum.tif"), slope = paste0(wd$data_p,
"temp_rasts/Slope.tif"), output = paste0(wd$data_p, "temp_rasts/SPI.tif"),
exponent = 1)
SPI <- rast(paste0(wd$data_p, "temp_rasts/SPI.tif"))
SPI <- project(SPI, DEM)
writeRaster(SPI, paste0(wd$data_p, "temp_rasts/SPI.tif"), overwrite = TRUE)
# SPI <- SPI * SaltLakeBinaryInvNA
SPI.01 <- RescaleSpatRast(SPI)
SPI.rs <- RescaleSpatRast(SPI.01) + 0.5
SPI.rs <- SPI.rs * WetAgInv
######### FreshLakeDepth #########
FreshLakeDepth[is.nan(FreshLakeDepth)] <- NA
FreshLakeDepth.12 <- RescaleSpatRast(FreshLakeDepth)
FreshLakeDepth.12 <- FreshLakeDepth.12 + 1
FreshLakeDepth.12[is.na(FreshLakeDepth.12)] <- 0
######### IntnsCost #########
IntnsCost = TRI.01 + SPI.rs + FreshLakeDepth.12
TRI.01.wetMod <- LakesBinaryInvNA * TRI.01
SPI.01.wetMod <- LakesBinaryInvNA * SPI.rs
IntnsCostNW = TRI.01.wetMod + SPI.01.wetMod
stz = st_as_stars(IntnsCost)
plt = ggplot() + geom_stars(data = stz) + scale_fill_viridis_c(option = "magma",
na.value = NA, name = "IntnsCost") + geom_sf(data = lake.s, fill = "turquoise1",
color = "turquoise1", alpha = 0.1) + geom_sf(data = CatchLims.sf, color = "white",
size = 0.5, alpha = 0) + ggtitle("Agricultural Intensification Cost Index in the SBOM",
"as a function of TRI, SPI and Freshwater Lake Depth") + coord_sf(datum = sf::st_crs(26914)) +
theme_bw() + labs(x = "Easting", y = "Northing")
ggsave("IntnsCostRast.png", plot = plt, device = "png", path = wd$figs, scale = 1,
width = 6, height = 5, units = "in", dpi = 1000)
plt
2.0.6 Agricultural Potential
Here we calculate a proxy for carrying capacity. To do so, we first rescale Irrigation Potential [1,2], Wetland AG potential [=2] and NPP [0.5,1.5] to assumed relative values for the BOM. Then we take the maximum raster cell value among all three, which serves as our raster of Agricultural potential.
The metrics calculated here are:
- AGPot.tot = Total agricultural potential = sum(max(WetAgPot, IrrigPot(TWI.), NPP))
- AGPot.avg = Mean agricultural potential = sum(max(WetAgPot, IrrigPot, NPP))/Area
- AGPotNW.tot = same as above but with 0 for wetland areas
- AGPotNW.avg = same as above but with 0 for wetland AG
######### AGPot #########
TWI.12.wetMod <- LakesBinaryInvNA * TWI.12
NPP.rs.wetMod <- LakesBinaryInvNA * NPP.rs
AGPot = max(WetAgPot2, TWI.12, NPP.rs)
AGPotNW = max(TWI.12.wetMod, NPP.rs.wetMod)
stz = st_as_stars(AGPot)
plt = ggplot() + geom_stars(data = stz) + scale_fill_viridis_c(na.value = NA,
name = "AGPot") + geom_sf(data = lake.s, fill = "turquoise1", color = "turquoise1",
alpha = 0.1) + geom_sf(data = CatchLims.sf, color = "white", size = 0.5,
alpha = 0) + ggtitle("Agricultural Potential Index in the SBOM", "As a function of NPP, TWI and Freshwater Lake") +
coord_sf(datum = sf::st_crs(26914)) + theme_bw() + labs(x = "Easting", y = "Northing")
ggsave("AGPotRast.png", plot = plt, device = "png", path = wd$figs, scale = 1,
width = 6, height = 5, units = "in", dpi = 1000)
plt
2.0.7 Population Pressure
With our proxy for carrying capacity in hand, we can now calculate a proxy for population pressure.
The metrics calculated here are:
- PopPressure = Population / Total Ag Potential = larger values indicate greater levels of pop pressure
- PopPressureNW = Pop Pressure with no wetland AG = Population / Total Ag Potential PopPressure
2.0.8 Erosion Potential
Landscape degradation in the form of erosion is thought to have been a major factor in settlement persistence across the premodern world. To create a proxy metric for erosion potential, we can use established laws of geomorphology and standard topographical proxies of their magnitude.
In quantitative geomorphology, erosion is classed into two major processes: Runoff Processes and Hillslope Diffusion Processes. Runoff is driven by overland flow while diffusion is long term soil creep or mass wasting.
The topographic variables diagnostic of runoff processes on hillslopes are slope hillslope length and catchment area. These are the topographic state variables in hillslope transport processes via continuous flowing water, such as sheetwash, rilling, and gullying. Runoff hillslope processes are dominant on middle and lower hillslopes, where water from uphill concentrates into progressively larger flows going downhill (Anderson & Anderson, 2010). Here, the sediment flux (or ‘discharge’), \(Q_{X}\), the amount of sediment transported through a point on a hillslope, can be represented by the equation \[ Q_{x}=kA^{m}\frac{\partial Z^{n}}{\partial X}=kA^{m}S^{n} \] where \(X\) is the point on the hillslope, \(Z\) is the depth of erosion, \(A\) is the upslope catchment area flowing through this point (assuming uniform rainfall and infiltration), \(S\) is the local slope, \(m\) and \(n\) are constants setting the relative importance of flowing water, and \(k\) is a parameter representing the erodibility of the sediment due to soil properties and surface conditions (Anderson & Anderson, 2010). Because \(m\) and \(n\) are >1 for processes dominated by concentrated runoff, this equation means that runoff and sediment discharge increase exponentially moving downhill. In three dimensions, \(A\) becomes the upslope catchment area flowing through a point in the hillslope. The rate of continuous runoff-based erosion at any point on the hillslope is given by the equation \[ \frac{\partial Z}{\partial t}=-\frac{1}{b} \frac{\partial Q_{X}}{\partial A} \] where \(b\) is the bulk density of the sediment and \(t\) is time (Anderson & Anderson, 2010). One proxy metric for topographic controls on hillslope runoff discharge and erosion dynamics is the Sediment Transport Index (STI), which provides a unit-less topographic index for potential erosion and deposition (Moore et al., 1993; Wilson & Burrough, 1999). Here, we calculate STI in R using the “whitebox” package (Lindsay, 2016).
The topographic variables diagnostic of hillslope sediment transport dynamics via diffusive processes are slope and curvature, the derivative of slope. These diffusive processes include rainsplash, bioturbation, tree-throw, discontinuous surface runoff, and mass wasting (from soil creep to landslides). Diffusion dynamics will be dominant wherever concentrated runoff is absent, so that diffusive sediment discharge is linearly dependent on slope alone. Slope is also a crucial threshold for hillslope stability in mass wasting events like landslides, which is frequently equated with the landscape’s observed angle of repose (Anderson & Anderson, 2010; Roering et al., 1999; Selby, 1993). The rate of erosion in hillslope diffusion processes is topographically controlled by curvature (the derivative of slope), such that \[ \frac{\partial Z}{\partial t}=-\frac{k}{b} \frac{\partial^{2} Z}{\partial X^{2}} \] where again \(k\) is the coefficient of diffusivity due to soil properties. Detailed empirical analyses have found that hilltop curvature is strongly correlated with erosion rates across transient landscapes (Hurst et al., 2012; Montgomery & Brandon, 2002; Palumbo et al., 2010). As such, the curvature of hillslopes and hilltops can provide crucial information on both local and landscape-wide erosion rates. Here, we calculate curvature in R using the “whitebox” package (Lindsay, 2016).
Given that runoff leads to much faster erosion than diffusion, we can scale Runoff Erosion Potential= STI to [1,2], and diffusion to [0,1.5], such that wherever runoff is prevalent it will be the dominant process. The max of these two rasters is then the preliminary erosion rate.
The rate of sediment and water transport through a landscape – and thus landscape degradation – are greatly impacted by its drainage density. Drainage density, \(D\), is the total length of a channel network, \(L\), divided by the area of its corresponding drainage basin, \(A\), such that for every stream \(i\) \[ D=\frac{\sum_{i = 1}^{n} L_{i}}{A} \] The denser the drainage network, the shorter hillslopes will be. Because the transport of both water and sediment are much faster and more efficient in channels than on hillslopes, higher drainage density results in higher erosion rates and rapid runoff discharge responses to rainfall. The drainage density of any given landscape will fluctuate based on the volume of runoff, which is primarily impacted by changes in vegetation density and rainfall intensity—especially in tropical environments. Drainage density is therefore a fundamentally empirical question, controlled by runoff discharge on hillslopes (Anderson & Anderson, 2010; Horton, 1945; Montgomery & Dietrich, 1988, 1992; Tucker & Bras, 1998; Tucker & Slingerland, 1997).
Nevertheless, topographic layout of different watersheds influences the potential density of their drainage networks. We can calculate a proxy for drainage density by using a simple flow accumulation raster in R calculated using the “whitebox” package (Lindsay, 2016). The max raster accumulation value is set at a level where we know streams are present (e.g. 10,000 = accumulation = upstream raster cells in contributing area). Then these values are rescaled [0,1], so that lower levels still contribute to the calculation but with less weight. Instead of drainage basins, we can simply sum the flow accumulation raster for each settlement catchment area (not the stream catchment, i.e. basin/watershed) and divide by the total settlement catchment area (with lake areas set at zero). This settlement catchment area drainage density value will be less than 1, making it the ideal denominator to multiply our preliminary erosion rate from above.
Given that Runoff Erosion Potential= STI[1,2] and Diffusion Erosion Potential= ProfCurv[0,1], the metrics calculated here are:
- ErosionPot.tot = Total Erosion Potential = sum(max(STI[1,2],Curv[0,1.5]))/DrainDens = take the max of the rescaled STI and Curv rasters, sum them over the catchment, and divide by the drainage density.
- ErosionPot.avg = Mean Erosion Potential = sum(max(STI[1,2],Curv[0,1]))/DrainDens/Area
######### STI #########
wbt_fd8_flow_accumulation(paste0(wd$data_p, "temp_rasts/DEM_breach.tif"), output = paste0(wd$data_p,
"temp_rasts/Accum_fd8.tif"), out_type = "specific contributing area")
Accum_fd8 <- rast(paste0(wd$data_p, "temp_rasts/Accum_fd8.tif"))
Accum_fd8 <- project(Accum_fd8, DEM)
writeRaster(Accum_fd8, paste0(wd$data_p, "temp_rasts/Accum_fd8.tif"), overwrite = TRUE)
wbt_sediment_transport_index(sca = paste0(wd$data_p, "temp_rasts/Accum_fd8.tif"),
slope = paste0(wd$data_p, "temp_rasts/Slope.tif"), output = paste0(wd$data_p,
"temp_rasts/STI.tif"))
STI <- rast(paste0(wd$data_p, "temp_rasts/STI.tif"))
STI <- project(STI, DEM)
writeRaster(STI, paste0(wd$data_p, "temp_rasts/STI.tif"), overwrite = TRUE)
STI = LakesBinaryInvNA * STI
STI.12 <- RescaleSpatRast(STI)
STI.12 <- STI.12 + 1
######### Curv #########
wbt_profile_curvature(dem = paste0(wd$data_p, "temp_rasts/DEM_breach.tif"),
output = paste0(wd$data_p, "temp_rasts/Curv.tif"))
Curv <- rast(paste0(wd$data_p, "temp_rasts/Curv.tif"))
Curv <- project(Curv, DEM)
writeRaster(Curv, paste0(wd$data_p, "temp_rasts/STI.tif"), overwrite = TRUE)
Curv = LakesBinaryInvNA * Curv
Curv[Curv < -0.001] <- -0.001
Curv.rs <- RescaleSpatRast(Curv) * 1.5
######### DrainDens Inputs #########
Accum2 <- Accum
Accum2 = LakesBinaryInvNA * Accum2
Accum2[Accum2 > 10000] <- 10000
Accum2 <- RescaleSpatRast(Accum2)
######### ErosionPot #########
ErosionPot <- max(STI.12, Curv.rs)
stz = st_as_stars(ErosionPot)
plt = ggplot() + geom_stars(data = stz) + scale_fill_viridis_c(option = "turbo",
na.value = NA, name = "ErosionPot") + geom_sf(data = lake.s, fill = "turquoise1",
color = "turquoise1", alpha = 0.3) + geom_sf(data = CatchLims.sf, color = "white",
size = 1, alpha = 0) + ggtitle("Topographic Erosion Potential Index in the SBOM",
"as a function of STI, Curvature and Drainage Density") + coord_sf(datum = sf::st_crs(26914)) +
theme_bw() + labs(x = "Easting", y = "Northing")
ggsave("ErosionPotRast.png", plot = plt, device = "png", path = wd$figs, scale = 1,
width = 6, height = 5, units = "in", dpi = 1000)
plt
rm(EZm, m, npp, stz, plt)2.1 Calculate Metrics for Catchments/Periods
Catch_List2 <- list() #create output list
for (i in 1:length(Catch_List)) {
# define catchment area polys as temp object
tmp.p <- Catch_List[[i]]
# convert to SpatVectors class for use with terra package (much much
# faster)
tmp.p2 <- vect(tmp.p)
# NPP
x <- terra::extract(NPP, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$NPP.tot <- x[, 2]
tmp.p@data$NPP.avg <- tmp.p@data$NPP.tot/tmp.p@data$Catchment_ha
# Env Zone
x <- terra::extract(EZ, tmp.p2, fun = mean, na.rm = T)
tmp.p@data$EZ.avg <- x[, 2]
x <- terra::extract(EZ, tmp.p2, fun = sd, na.rm = T)
tmp.p@data$EZ.sd <- x[, 2]
# TRI (Terrace Need)
x <- terra::extract(TRI, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$TRI.tot <- x[, 2]
tmp.p@data$TRI.avg <- tmp.p@data$TRI.tot/tmp.p@data$Catchment_ha
## Hydraulic Agriculture Potential Irrigation Potential
x <- terra::extract(TWI.01, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$IrrigPot.tot <- x[, 2]
tmp.p@data$IrrigPot.avg <- tmp.p@data$IrrigPot.tot/tmp.p@data$Catchment_ha
# Wetland Agriculture Potential
x <- terra::extract(WetAgPot, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$WetAgPot.tot <- x[, 2]
tmp.p@data$WetAgPot.avg <- tmp.p@data$WetAgPot.tot/tmp.p@data$Catchment_ha
# Intensification costs
x <- terra::extract(IntnsCost, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$IntnsCost.tot <- x[, 2]
tmp.p@data$IntnsCost.avg <- tmp.p@data$IntnsCost.tot/tmp.p@data$Catchment_ha
# Intensification costs excluding wetland agriculture
x <- terra::extract(IntnsCostNW, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$IntnsCostNW.tot <- x[, 2]
tmp.p@data$IntnsCostNW.avg <- tmp.p@data$IntnsCostNW.tot/tmp.p@data$Catchment_ha
# Agricultural Potential
x <- terra::extract(AGPot, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$AGPot.tot <- x[, 2]
tmp.p@data$AGPot.avg <- tmp.p@data$AGPot.tot/tmp.p@data$Catchment_ha
# Agricultural Potential excluding wetland agriculture
x <- terra::extract(AGPotNW, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$AGPotNW.tot <- x[, 2]
tmp.p@data$AGPotNW.avg <- tmp.p@data$AGPotNW.tot/tmp.p@data$Catchment_ha
# Population Pressure
tmp.p@data$PopPressure <- tmp.p@data$Population.s2/tmp.p@data$AGPot.tot
tmp.p@data$PopPressureNW <- tmp.p@data$Population.s2/tmp.p@data$AGPotNW.tot
# Erosion Potential
x <- terra::extract(ErosionPot, tmp.p2, fun = sum, na.rm = T)
tmp.p@data$ErosionPot.tot <- x[, 2]
tmp.p@data$ErosionPot.avg <- tmp.p@data$ErosionPot.tot/tmp.p@data$Catchment_ha
Catch_List2[[i]] <- tmp.p #save to output list
}
# rename list items
catch.names <- names(Catch_List)
names(Catch_List2) <- catch.names
Catch_List <- Catch_List2
rm(Catch_List2)2.2 Reorganize and Export Raster Data
TopoEnvData_RastStack <- c(DEM, Curv, Curv.rs, Accum, Accum_fd8, Accum2, AGPot,
AGPotNW, EZ, FreshLakeDepth, FreshLakeDepth.12, IntnsCost, IntnsCostNW,
LakesBinaryInv, LakesBinaryInvNA, NPP, NPP.rs, NPP.rs.wetMod, SaltLakeBinaryInv,
SaltLakeBinaryInvNA, slope, SPI, SPI.01, SPI.01.wetMod, SPI.rs, STI, STI.12,
TRI, TRI.01, TRI.01.wetMod, TWI, TWI.01, TWI.12, TWI.12.wetMod, WetAgInv,
WetAgPot, WetAgPot2, ErosionPot)
nam = c("DEM", "Curv", "Curv.rs", "Accum", "Accum_fd8", "Accum2", "AGPot", "AGPotNW",
"EZ", "FreshLakeDepth", "FreshLakeDepth.12", "IntnsCost", "IntnsCostNW",
"LakesBinaryInv", "LakesBinaryInvNA", "NPP", "NPP.rs", "NPP.rs.wetMod",
"SaltLakeBinaryInv", "SaltLakeBinaryInvNA", "slope", "SPI", "SPI.01", "SPI.01.wetMod",
"SPI.rs", "STI", "STI.12", "TRI", "TRI.01", "TRI.01.wetMod", "TWI", "TWI.01",
"TWI.12", "TWI.12.wetMod", "WetAgInv", "WetAgPot", "WetAgPot2", "ErosionPot")
names(TopoEnvData_RastStack) = nam
terra::writeRaster(TopoEnvData_RastStack, paste0(wd$data_p, "TopoEnvData_RastLayers.tif"),
overwrite = T)
terra::writeRaster(TopoEnvData_RastStack, paste0(wd$data_f, "TopoEnvData_RastLayers.tif"),
overwrite = T)
unlink(paste0(wd$data_p, "temp_rasts"), recursive = TRUE)
rm(Catch_List2, x, DEM, Curv.01, Curv, Curv.rs, Accum, Accum_fd8, Accum2, AGPot,
AGPotNW, EZ, FreshLakeDepth, FreshLakeDepth.12, IntensCost, IntnsCostNW,
LakesBinaryInv, LakesBinaryInvNA, NPP, NPP.01, NPP.rs, NPP.rs.wetMod, NPP2,
SaltLakeBinaryInv, SaltLakeBinaryInvNA, slope, SPI, SPI.01, SPI.01.wetMod,
SPI.rs, STI, STI.12, tmp.p, tmp.p2, tmp.r, TRI, TRI.01, TRI.01.wetMod, TRI.12,
TRI.rs, TWI, TWI.01, TWI.12, TWI.12.wetMod, WetAgInv, WetAgPot, WetAgPot2,
ErosionPot, IntnsCost, npp, EZm, TopoEnvData_RastStack)3 Recombining and Reorganizing the Data
# Convert lists of period-wise sites/catchments to single SPDF objects
All_Agg_SitePoly <- do.call(rbind, Poly_List)
All_Agg_CatchPoly <- do.call(rbind, Catch_List)
# variables from site data that needs transfer over to catchment areas
###colz1 = setdiff(colnames(All_Agg_SitePoly@data),colnames(All_Agg_CatchPoly@data))
# variables from catchment areas that needs transfer over to site data
colz2 = setdiff(colnames(All_Agg_CatchPoly@data),colnames(All_Agg_SitePoly@data))
#reorder the data to match
All_Agg_SitePoly <- All_Agg_SitePoly[order(All_Agg_SitePoly$AggSite),]
All_Agg_CatchPoly <- All_Agg_CatchPoly[order(All_Agg_CatchPoly$AggSite),]
#check to see that the two datasets are in the right order
#identical(All_Agg_SitePoly@data$AggSite, All_Agg_CatchPoly@data$AggSite)
# Site data to catchment areas
###Site_to_Catch <- All_Agg_SitePoly@data %>% dplyr::select(!!!syms(colz1))
###All_Agg_CatchPoly@data <- cbind(All_Agg_CatchPoly@data,Site_to_Catch)
# catchment area data to sites
Catch_to_Site <- All_Agg_CatchPoly@data %>% dplyr::select(!!!syms(colz2))
All_Agg_SitePoly@data <- cbind(All_Agg_SitePoly@data,Catch_to_Site)
#Reorganize the data
ordering <- c(
#ID VARIABLES
"AggSite","AggID","Site","East","North","SurvReg","Number","CerPhase","Period",
"PeriodType","PeriodLength","PeriodNum","PeriodBegin","PeriodEnd",
"OccSeqLoc","OccSeqLoc.Sites","SubOccSeqLoc","SubOccSeqLoc.Sites",
"ComponentNum", "ComponentSites",
#CHRONOLOGICAL VARIABLES
"PeriodInterval", "PeriodBegin", "PeriodBegin.era", "PeriodMidpoint",
"PeriodMidpoint.era", "PeriodEnd", "PeriodEnd.era", "PeriodLength",
#OCCUPATION VARIABLES (COUNTS)
"Occ.EF","Occ.EF_MF","Occ.MF","Occ.MF_LF","Occ.LF","Occ.LF_TF",
"Occ.TF","Occ.TF_CL","Occ.CL","Occ.CL_ET","Occ.ET","Occ.ET_LTAzI",
"Occ.LTAzI","Occ.LTAzI_EA","Occ.EA","Occ.EA_LA","Occ.LA","Occ.TOT",
#SUBOCCUPATION VARIABLES (COUNTS)
"SubOcc.EF","SubOcc.EF_MF","SubOcc.MF","SubOcc.MF_LF","SubOcc.LF",
"SubOcc.LF_TF","SubOcc.TF","SubOcc.TF_CL","SubOcc.CL","SubOcc.CL_ET",
"SubOcc.ET","SubOcc.ET_LTAzI","SubOcc.LTAzI","SubOcc.LTAzI_EA",
"SubOcc.EA","SubOcc.EA_LA","SubOcc.LA","SubOcc.TOT",
#SITE AREA AND OCCUPATIONAL DENSITY VARS
"Area_ha","Perim_m2","SherdDens","Tot.Assemb","FwOvlp.Assemb",
"BwOvlp.Assemb", "Net.Assemb",
#STEP #2 DEMOGRAPHIC VARIABLES
"Population.s2","Log_Population.s2", "ApportAssemb", "PopDens.s2",
"UrbanScale.s2", "UrbanPop.s2","RuralPop.s2", "PctUrban.s2","PctRural.s2",
"Prior", "Observed", "MeanOccuProb",
#STEP #2 DEMOGRAPHIC RATES
"Pct_deltaPop12", "Pct_deltaPop01", "r12_Pert", "r01_Pert",
"r23_Pert", "rPert", "rPert0", "rPert2",
"Pct_UrbdeltaPop12", "Pct_UrbdeltaPop01", "Urb_r12_Pert", "Urb_r01_Pert",
"Urb_r23_Pert", "Urb_rPert", "Urb_rPert0", "Urb_rPert2",
"Abs_rPert","RS_rPert","LogRS_rPert","LogAbs_rPert","Q.rPert","Q.Abs_rPert",
"Q.RS_rPert","Q.LogRS_rPert","Q.LogAbs_rPert","Urb_Abs_rPert",
"Urb_RS_rPert","Urb_LogRS_rPert","Urb_LogAbs_rPert","Q.Urb_rPert","Q.Urb_Abs_rPert",
"Q.Urb_RS_rPert","Q.Urb_LogRS_rPert","Q.Urb_LogAbs_rPert",
"Qp.rPert","Qp.Abs_rPert","Qp.RS_rPert","Qp.LogRS_rPert","Qp.LogAbs_rPert",
#STEP #1 DEMOGRAPHIC VARIABLES
"Population.s1","PopDens.s1","UrbanScale.s1","UrbanPop.s1","RuralPop.s1",
"PctUrban.s1","PctRural.s1",
#CONTINUITY VARIABLES
"AreaBwCont","AreaFwCont","PopBwCont","PopFwCont","FwOvlp.Sites",
"FwOvlp.Area","FwOvlp.Pop","BwOvlp.Sites","BwOvlp.Area","BwOvlp.Pop",
#PERSISTENCE VARIABLES
"Found","FoundInit","Abandon","Persist","DewarType",
"OccuTime","OccuInertia","UrbOccuTime","UrbOccuInertia",
#STEP #4 TRANSPORT NETWORK VARIABLES
"TranspDens.pct", "TranspDens.rank","centr_deg","centr_btw","centr_eig",
"centr_clos","centr_hrmo","centr_hub","centr_auth","centr_pgrk",
"centr_deg_n","centr_btw_n","centr_eig_n","centr_clos_n","centr_hrmo_n",
"centr_hub_n","centr_auth_n","centr_avg","trans_loc","trans_locavg",
"density","connectiv","trans_glob","cntrlz_deg","cntrlz_btw","cntrlz_eig",
"cntrlz_clo","cntrlz_hub","cntrlz_auth","cntrlz_deg_n","cntrlz_btw_n",
"cntrlz_eig_n","cntrlz_clo_n","cntrlz_hub_n","cntrlz_auth_n","cntrlz_avg",
#STEP #4 CATCHMENT AREA AND POP DENSITY VARIABLES
"Catchment_ha", "CatchmentBeyond_ha", "Catch_Popdens",
"Catch_Popdens_PropMax","Catch_Popdens_Rank","CatchB_Popdens",
"CatchB_Popdens_PropMax", "CatchB_Popdens_Rank",
#STEP #5 CATCHMENT ENVIRONMENT/TOPOGRAPHY VARIABLES
"NPP.tot","NPP.avg","EZ.avg","EZ.sd","TRI.tot","TRI.avg","IrrigPot.tot",
"IrrigPot.avg","WetAgPot.tot","WetAgPot.avg","IntnsCost.tot",
"IntnsCost.avg","IntnsCostNW.tot","IntnsCostNW.avg","AGPot.tot",
"AGPot.avg","AGPotNW.tot","AGPotNW.avg","PopPressure","PopPressureNW",
"ErosionPot.tot","ErosionPot.avg",
#SURVEY METADATA
"M_Sites","M_SiteCode","M_SiteName","M_FieldSite.Region",
"M_FieldSite.Period","M_SurveyYearNumber","M_Supervisor","M_Map",
#OLD tDAR BOM SURVEY VARIABLES
"O_Elev","O_ElevMed","O_ElevMin","O_ElevMax","O_EZcode",
"O_EnvironmentalZone","O_Soil","O_SoilMed","O_SoilMin","O_SoilMax",
"O_Erosion","O_ErosionMed","O_ErosionMin","O_ErosionMax","O_ModernUse",
"O_ModernSettlement","O_Rainfall","O_Area","O_MoundDomestic",
"O_MoundCeremonial","O_MoundQuestionable","O_MoundTotal",
"O_MoundRecorded","O_DMoundArea","O_Architecture","O_TerraceConfidence",
"O_TerraceExtent","O_Sherd","O_SherdMed","O_SherdMin","O_SherdMax",
"O_Rubble","O_RubbleMed","O_RubbleMin","O_RubbleMax","O_Population",
"O_PopMin","O_PopMax","O_PopMethod","O_stcode","O_SiteType",
"O_SubPeriod1","O_SubPeriod2","O_OccEF","O_OccMF","O_OccLF","O_OccTF",
"O_OccCL","O_OccEC","O_OccMC","O_OccLC","O_OccET","O_OccLT","O_OccAZ",
"O_OccEA","O_OccLA","O_OccTot","O_OccSeqLoc","O_SubOc1","O_SubOc2",
"O_PdDupSite","O_Group","O_Comments")
#Make sure everything is kosher
#setdiff(colnames(All_Agg_CatchPoly@data),ordering)
#setdiff(colnames(All_Agg_SitePoly@data),ordering)
#setdiff(ordering,colnames(All_Agg_CatchPoly@data))
#setdiff(ordering,colnames(All_Agg_SitePoly@data))
# Reorder the data for both sites and catchment areas
All_Agg_SitePoly@data <- All_Agg_SitePoly@data %>% dplyr::select(!!!syms(ordering))
All_Agg_CatchPoly@data <- All_Agg_CatchPoly@data %>% dplyr::select(!!!syms(ordering))4 Export Data for Script #6
# AggSite polygons
writeOGR(All_Agg_SitePoly, paste0(wd$data_p, "SBOM_AggSitePoly5.gpkg"), "SBOM_AggSitePoly5",
driver = "GPKG", overwrite_layer = TRUE)
# Catchment areas
writeOGR(All_Agg_CatchPoly, paste0(wd$data_p, "SBOM_CatchPoly5.gpkg"), "SBOM_CatchPoly5",
driver = "GPKG", overwrite_layer = TRUE)